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We discuss the effects of external stress on the thermal denaturation of homogeneous DNA. Pulling 
double-stranded DNA at each end exerts a profound effect on the thermal denaturation, or melting, 
of a long segment of this molecule. We discuss the effects on this transition of a stretching force 
applied to opposite ends of one of the DNA strands, including full consideration of the consequences 
of excluded volume, the analysis of which is greatly simplified in this case. We also discuss the 
interplay of thermal denaturation and force-generated separation when the tension is generated by 
a force at the end of the duplexed strands and an equal and opposite force is applied to the other 
end of the second strand. 

PACS numbers; 



I. INTRODUCTION 

When heated, double-stranded DNA develops regions 
of strand separation, known as "denaturation bubbles." 
[l[ As the temperature increases, these bubbles grow in 
size and proliferate. The culmination of this process is 
complete separation of the strands, that is the complete 
thermal denaturation — or melting — of DNA, a transition 
of significant biological and technological significance. 
Among the most important theoretical approaches to this 
process is a collection based, or mathematically related, 
to an underlying model introduced by Poland and Scher- 
aga [3, Q . These approaches include the formalism intro- 
duced by Peyrard and Bishop j^] that maps the process of 
melting into the disappearance of a bound state in a one- 
dimensional potential well. Given the results obtained 
with the use of this family of models, the melting process 
in the case of homogeneous double-stranded DNA — that 
is, the molecule consisting of, say, one strand containing 
only cytosines and another containing only guanines — 
is well-understood. It is known with a relatively high 
degree of certainty that an infinitely long molecule of 
this kind will undergo a first order melting transition Q . 
This conclusion follows from the consideration of the ef- 
fects of self-avoidance — in particular, the consequences of 
self avoidance with regard to the structure of the vertex 
connecting an intact portion of DNA with a denatura- 
tion bubble. The treatment of this process follows from 
the work of Duplantier [1, 0] on the renormalization of 
vertices for arbitrary polymer networks. 

The sharp transition to the completely denatured state 
is an example of a phase transition in a one-dimensional 
system. The existence of such a transition follows from 
the effective long-range interaction inherent in the sta- 
tistical mechanics of the denaturation bubbles. The cor- 
relations that are propagated in a bubble result in an 
effective Boltzmann factor consistent with the statisti- 
cal mechanics of an inverse square Ising model, which is 
known to undergo a phase transition [1, Q . 

In light of the expectation of a sharp transition in ther- 
modynamic limit in the case of homogeneous DNA, it 
is noteworthy that experiments on biological DNA pro- 



duce melting curves that are belie the expectation of a 
true phase transition in the thermodynamic limit. Here, 
we utilize the standard definition of a phase transition 
as non-analyticity in thermodynamic functions. What 
is seen experimentally is, rather a collection of highly- 
structured, but nevertheless smooth, curves when for 
instance, the specific heat is measured (see, e.g. lu]). 
It appears that the inhomogeneity inherent in the DNA 
present in living organisms, which effectively translates 
to random inhomogeneity in the context of thermal de- 
naturation, profoundly affects the nature of the transi- 
tion. This is consistent with the fact that such random 
inhomogen eity is relevant in the sense of the Harris cri- 
terion [II], [13] ■ A brief demonstration that this is so is 
presented in Appendix [Cl 

An additional mechanism for the separation of DNA 
strands is the application of external forces, which under 
appropriate circumstances leads to the "unzipping" (l3| 
or to the stretching-induced "melting" of the molecule. 
This process has been investigated theoretically [13, El, 
[TgI . [TtI . [isj as has the interplay between the effects of ex- 
ternally applied force and thermal fluctuations. In the 
latter case, the melting process is modeled as a helix-coil 
transition, which can be represented in terms of a sys- 
tem with an intimate relationship to the one-dimensional 
Ising model with short range interactions. In this ap- 
proach to the denaturation process, there is no prospect 
of a thermodynamic transition p^ . 

The set of calculations that we report here assumes 
the underlying validity of the Poland-Scheraga-based ap- 
proaches to DNA melting. We consider melting as the 
result of the accumulation and possible merging of de- 
naturation bubbles. What is added to the picture is a 
pair of equal and opposite forces at the two ends of the 
denaturing strands. Restricting our focus to homoge- 
neous DNA, we assess the consequences of this force pair 
under two circumstances: 

a: The forces are applied either to both strands simul- 
taneously or to both ends of one of them. 

b: The force at one end acts on one strand, while the 
force at the opposite end acts on the other strand. 
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These two ways of pulling at the ends of duplexed DNA 
are illustrated in Fig. [T] 




b 

FIG. 1: The two ways in which a strand of duplexed DNA can 
be stretched, either a: with forces applied to opposite ends 
of the same strand, or b: with a force at the left end of one 
of the strands and an equal and opposite force applied to the 
right end of the other strand. 



energy cost of a self intersection in light of the forces act- 
ing on the ends of the chain. In Section IIIII we focus on 
the one way in which self- avoidance influences the statis- 
tical mechanics of the melting transition, in the interior 
of a denaturation bubble. We find that it introduces 
a logarithmic modification to the mean field result for 
melting exponents. Section IIVI is devoted to a general 
analysis of the infiuence of self-avoidance and stress on 
thermal denaturation of DNA, particularly as it relates 
to the phase diagram and key temperature dependences. 
In Section|V]we present a preliminary analysis of the con- 
sequences of stress applied as in Fig. [T}). We find that 
pulling on the separate strands renders the undenatured 
state metastable at any temperature. If there is a tran- 
sition, it is of the spinodal type [2l'| . Finally, in Section 
I VII we review our analysis and point to its implications. 



In the first case above, the effect of the force is two- 
fold. First, as noted by Rouzina and Bloomfield pTf . 
stretching both strands produces an energetic bias in fa- 
vor of denaturation, in that single strands of DNA are 
more easily elongated than the intertwined strands of the 
duplexed version of the molecule. On the other hand, 
DNA melting represents the classic competition between 
energy and entropy, with denaturation bubbles embody- 
ing the entropically favored, energetically costly state. 
Stretching a pair of strands has the effect of applying an 
additional energy penalty to the states available to the 
separated strands. Thus, the imposition tension on melt- 
ing strands of DNA can either promote or inhibit thermal 
denaturation. As we will see, both consequences can be 
observed. 

The second version of an externally generated tension 
has the unequivocal effect of favoring strand separation. 
In fact, it is straightforward to demonstrate that pulling 
one strand at one end and the other strand at the op- 
posite end causes the state of the molecule in which any 
portion is duplexed to be a metastable state, rather than 
a state of true equilibrium. 

The paper proceeds as follows. In Section |ll] we re- 
view the effect of stress on the statistics of a Gaussian 
chain of equal and opposite forces applied to both ends 
and introduce the mathematical tools for the analyzing 
the phase transition in terms of the nonanalyticity of the 
generating functions. In particular, we recall that the 
asymptotic statistics of a system with a fixed number 
of monomeric units is controlled in the thermodynamic 
limit by the singularity in the generating function lying 
closest to the origin in the complex fugacity plane [20|. 
Given this, we are able to show that almost all modi- 
fications of the generating function that follow from ex- 
cluded volume considerations will exert a negligible effect 
on the melting transition, in that the singularities asso- 
ciated with those modifications are further away from 
the origin than the singularity that arises in the case of 
the unrestricted chain. The mitigation of these excluded 
volume effects can be simply understood in terms of the 



II. THE EFFECT OF AN EXTERNAL FORCE 
ON A FREELY-JOINTED CHAIN 

For the purposes of calculating the effect of the ex- 
ternal forces on a single strand of DNA, an appropriate 
starting point is the freely jointed chain (FJC). A long 
segment of DNA can be approximated as a chain of many 
identical molecules connected with each other at joints 
which allow for spatial rotations (See Fig. [2]). 




FIG. 2: The Gaussian chain in which the monomeric units are 
separated by a distance a can be modeled as a freely jointed 
chain with link length I. I \s the Kuhn length of the chain. 



The probability distribution of the end-to-end distance 
vector R of such an object is given by 
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where N is the number of units of the FJC, I is the length 
of each unit. It can be shown that I is equal to twice the 
persistence length of the chain . 

The first delta function on the right hand side of Eq. ([T|) 
ensures that the vectors Z„ of the chain elements add up 
to the distance vector R. 

If a force F is applied at one end of the chain and 
an equal and opposite force acts on the other end, then 
there is an additional Boltzmann factor weighting the 
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chain configuration sum equal to exp[F ■ R/ksT] so that 
Eq. Il]) becomes 



Pn{R,F) = exp 
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Using the saddle point approximation, the distribution 
function in Eq. ([2]) can be evaluated as 
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where we have introduced the dimensionless parameters 
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/ is the magnitude of /; p\\ and pj_ are the components 
of p correspondingly parallel and perpendicular to the 
dimensionless force /, and 



1 



p\ = coth/ - -. 



(5) 



The details of the calculation are presented in Ap- 
pendix]^ We have introduced the definitions of the func- 
tions '^flpAf) and «'^/p||(/) in Eqs. jUni and dMH). 
In the limit of weak forces, / ^ 1, 
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and the distribution function in Eq. Q splits into a prod- 
uct of two terms, one describing the Gaussian chain and 
the other the statistical weight due to a small-force per- 
turbation: 
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To evaluate the end-to-end separation of the chain under 
the influence of the external force we compute (i?^), 



J{Rj + RI)Pn{R, f)dR^^dRj 
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In the limiting case of weak forces / ^ 1 from Eq. 
with the use of Eq. it follows 
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Thus the DNA chain behaves as a Gaussian coil, with 
r.m.s. end-to-end distance proportional to ^/N. 
For strong forces / S> 1, 
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that corresponds to a fully stretched chain. 

We acquire more essential information by constructing 
the Fourier transformed generating function G(z,q,/), 
defined formally in terms of a sum over monomeric units 
k of the spatial Fourier transform of configurations. 

We find 



G{z,q,f) = Y,^'Pk{q,f) 
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where 
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is the critical fugacity, with 

=ap|>0. 
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and a is the size of a monomeric unit. 

The last line of Eq. (fTTj) reflects the fact that we are in- 
terested in the behavior of the generating function when 
z zq and in the limit of small q. 

We will extract from the above generating function 
the number of weighted configurations of a /c-monomer 
chain. We do this by calculating the coefficient of z'^ 
in the power series expansion of the generating function, 
which we take to be given by the last line of Eq. (|lip . 
The actual calculation makes use of Cauchy's theorem 
[23 II and involves the following contour integration: 
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where the contour is as illustrated on the upper left hand 
side of Fig. [3l The evaluation of the integral involves the 
distortion of the integration contour, as indicated in the 
lower right hand side of Fig. [3l so as to enclose the singu- 
larities in the generating function. The figure illustrates 
the result of that distortion when the singularities are 
two simple poles. In our case, the generating function 
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FIG. 3: The contours utilized in the evaluation of the integral 
Eq. p6|) leading to the extraction of the k^^ power of z in the 
expansion of the generating function. On the upper left hand 
side: the original contour, consistent with the extraction of 
that coefficient via Cauchy's theorem. On the lower right 
hand side: the distortion of the original contour to enclose 
singularities of the integrand — in the illustrated case simple 
poles — located at the heavy dots in the figure. 



possesses a single pole, which lies on the real axis. When 
there is more than one singularity, the dominant contri- 
bution arises from the singularity that lies closest to the 
origin. In fact, in the thermodynamic limit, k ~* oo, 
effectively the only contribution that matters is the one 
generated by the singularity closest to 2; = 0. This is not 
an issue in the calculation performed here, but it will be 
as we consider the mathematics of melting as embodied 
in the Poland Scheraga model and modifications thereof. 

Continuing, we note, as indicated immediately above, 
that the contour integration is dominated by the single 
contribution at the pole in the function on the last line 
of Eq. (|16p. with q = 0, corresponding to the solution of 
the equation 



l-z/zo = 0. 



(17) 



Inserting this solution into the result of the contour in- 
tegration illustrated in the lower right hand corner of 
Fig. [3l and making note of the residue, we end up with 
the following result for the total number of weighted k- 
monomer configurations when the polymer is subjected 
to the externally generated tension / 



1 



G{z,q^O,f] 



dz 



/sinh/^ 



(18) 



Corrections for excluded volume: one loop 
order 



We can now assess the interplay of this externally- 
generated tension and self-avoidance in influencing the 
asymptotic statistics of an excluded volume Gaussian 
polymer. A way to assess this interplay is to consider 
the lowest order correction to the generating function, as 
shown in Fig. H] 




FIG. 4: The first order correction to the generating function 
arising from excluded volume. The dashed line is the effective 
repulsive interaction arising from self-avoidance. 



The one-loop correction corresponds to the expression 

u [ G{z,q,f )dq\\dqx_ 

1 ~ z/ zq + X\\q^^ + X^q\ - iY\\q\\ dq\\dq±_ 
d^q 
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where u is the coupling constant measuring the strength 
of the repulsive interaction. A{z) is a singularity-free 
function of z. The next-to-last line in Eq. ([TO]) follows 
from a shift in the contour of integration over the com- 
ponent of q parallel to /. 

The outcome of this calculation is that the one-loop 
correction to the effective self-energy of the generating 
function of the excluded volume yields an expression hav- 
ing the following form 



G(z,q,/) = 



1 - z/zo + X\\qf + X^ql - iY\iq\\ 



+ uA{z) - u^l - z/zo + y||^/4X|| 



(20) 



Setting q = 0, in order to locate the singularity that 
dominates a calculation of the total number of weighted 
configurations, we find 



G{z,q = 0,f) = 



1 — z/zq + u A(z) 
-1 
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The singularity in the function, if u is sufficiently small, 
is slightly shifted from the unrestricted value. 
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However, the singularity of the one-loop contribution 
is given by 



1 + 



In the weak-force limit it is equal to 

la 



(22) 



(23) 



In particular, the singularity in the generating function 
remains closer to the origin than the singularity in the 
one-loop contribution in Eq. This means that the 

latter singularity will not have an effect on the asymp- 
totic statistics of the polymer chain under tension, in 
contrast to the situation of the unstretched chain, for 
which the one-loop correction for self-avoidance exerts 
an essential, and transforming, influence on asymptotic 
configurational statistics. The shift is indicated in Fig.O 
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FIG. 5: At lowest nontrivial loop order, the original singu- 
larity in the generating function for q = 0, at « = zo, (black 
circle) is shifted a distance 0(u) to the right (cross), and lies 
closer to the origin than the singularity in the one-loop cor- 
rection at 2 = «c (open circle). 



B. Two loop order correction 

The next obvious question is whether this difference 
between the singularity structure of the corrected gen- 
erating function and the singularities of the corrections 
persists to higher loop order. A cursory investigation of 
the general structure of such contributions to asymptotic 
configurational statistics indicates that, if the strength 
of the self-avoiding interaction is low enough, the singu- 
larity structure of higher loop order corrections will not 
lead to an essential renormalization of loop statistics. 

As the argument utilized here is based on the structure 
of excluded volume corrections in real space, we begin 
with the form of the generating function in real space, 
the inverse Fourier transform of G'(z, ,q±,f). We find 
for this quantity 



G{z,R,f) 
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where we have introduced a new variable 
r — (r_L, z) = 



(25) 



In particular, in the limit of weak forces, / <C 1, we 
find in the frist order of / 



X|| = X_L = la/6, Y\\ = af/5, 



and 



G{z,r,f) cx 
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(26) 



(27) 



where we assumed that the force in the z direction. 

The real space version of the one-loop self energy 
shown in Fig.[4]can be reproduced to within non-singular 
contributions by allowing r to become very small in 
Eq. (|27p . looking in particular at the zeroth order in r 
terms in the resulting expression. Expanding the result 
in powers of r, we find for the one- loop correction 

- + vW^Vl - ^/^c + + O (r) , (28) 



where Zc is given by Eq. 

Note that the next to last term in Eq. (pS)) depends 
on the direction in which the displacement vector points. 
This contribution is, in fact, a singular one, in that it is 
indeterminate in the limit r — 0. However, the averag- 
ing inherent in compensating for self-intersection is over 
directions, which means that we must take the angular 
average {cos 6), which is equal to zero. We thus recover 
the singular portion of the one- loop correction as propor- 
tional to -y/l — z/zc, as we did in the previous section. 

We now turn to the two loop self energy. Fig. [6] illus- 
trates the most interesting two loop term. The structure 




FIG. 6: The two loop insertion. 



of this figure is fairly straightforward. It has the following 
spatial Fourier transform 



J e"^ ''G(z,r,/ = 0)3dV. 



(29) 
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As the dominant issue for us is the way in which this 
correction affects the total number of configurations, we 
are interested in the q = limit of the above expression. 
Inserting Eq. ([?f|) into Eq. with / set equal to zero, 
we obtain the following result for the two loop self energy 



~r^Jl-z/zo+pa/6l^6/al \ 3 



(30) 



When there is tension, the integrand in Eq. ([50]) is 
replaced by the one that contains a Boltzmann factor as- 
sociated with the applied forces. The integral to perform 
in this case is 



(X u~ I dfl / r dr 

Ira 
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(31) 



Carrying out the integrations, we end up with the fol- 
lowing result for the two loop self energy associated with 
the diagram in Fig. [Sj 



{g [(avWv'l - z/zo + pa/6l - f/l) ro 
i^e/^Wl - z/zo + Pa/6l + f/i) ro] } , (32) 



where 



g{x) ==e-^ + j;Ei(-x). 



(33) 



Taking the ro ^ limit of Eq. (|32|) , we find that there is 
a singularity of the form w In it;, where 



w = SVe/o^Vl - z/zQ + pa/m ± f/l. (34) 

This tells us that the leading order singularity in the two 
loop self energy lies on the following location on the real 
z axis: 



z = zo{l + 4 fa/271). 



(35) 



This is to the left of Zc = zo{l + pa/61), but further 
from the origin than zq, the leading singularity of the 
unrestricted linear chain polymer under the influence of 
tension, as given by Eq. (fT2|l . 

A simple — but we believe essentially correct — 
argument holds that the energy penalty associated with 
the necessity of looping back against the tension that is 
required for self-intersection militates against a funda- 
mental alteration of configurational statistics when cor- 
rections are made for excluded volume. For a more ex- 
tended discussion, see Appendix |BJ 



III. EFFECT OF SELF- AVOIDANCE ON THE 
DENATURATION BUBBLE 



of the system. This is through the modification of the 
contribution of the denaturation bubble to the partition 
function of the system. As noted previously, configura- 
tions in which self-intersection requires that the chain 
loop back on itself, in opposition to the imposed ten- 
sion, are, by reason of the energy penalty associated with 
such a configuration, rendered irrelevant to the asymp- 
totic statistics of thermal denaturation. This means that 
the vertex correction found by Kafri et al. Q to cause 
melting to become first order in the absence of tension 
plays no such transformative role when that tension is 
present. However, intersections of the two strands in the 
denaturation bubble as shown in Fig. [7| must be taken 
into account in the evaluation of the partition function 
sum of the melting chain. 



FIG. 7: The denaturation bubble, with corrections for self- 
interaction incorporated. Dotted lines show the repulsive in- 
teractions between the strands. 



In the denaturation bubble there are equal number of 
units in its lower and upper chains. To enforce this con- 
dition in the context of the calculation that we perform 
below, we assign different fugacites, zi^2 to the units in 
the two strands in the bubble. The relation between re- 
sulting generating function of the, in general, asymmetric 
loop ^^(zi, Z2)and the generating function representing 
the denatured loop is the following 

n(.) = ^/^^(«.Vi,^). 

Ztti J w w 

(36) 

where the contour of integration is around the origin. 
Indeed, by definition. 



.F(zi,Z2)=5ZEC'™n^"^2, 



(37) 



where Cnm is the partition function of the loop formed 
by two chains of n and m units correspondingly. 

Then according to Eq. , the generating function of 
the denaturation bubble becomes 



2/2 n/2 1 



t { m — n 1 
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y^ CmnZ"^^^ z"'^'^6rnn = C'rm^", (38) 



In the case of the melting DNA chain, there is one 
way in which self avoidance exerts a fundamental influ- 
ence on the statistics, and thereby the thermodynamics. 



that shows that n(z) is indeed the generating function of 
a bubble with equal number of units in each chain. In the 
second line of Eq. ([55)1 we have deformed the contour of 
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integration to a unit circle with the center at the origin. 
To account the possible interactions of the two chains we 
need to perform a Dyson-like summation over the loops 
formed by the strands interactions as shown in Fig|71 
Thus the generating function for the self-interacting de- 
naturation bubble becomes 

n(z) = :^/ — f:(-<^"(-V^,^). (39) 

l-Kl j W ^ — ' W 

n 

In the next section we derive the closed-form expression 
for n(z). 



The angular integral in Eq. (|4T|) leads to 



-air _ „-02r 



(43) 



where we have introduced the following notations 



+ y^l-Z2/io + 2/f|/4a;|| T y\\/V^- (44) 



Multiplying Eq. by and integrating from r = to 
oo we obtain 



oo -air 



-dr. 



(45) 



The integral in Eq. (|l5|) can be rewritten as 



A. Mutual avoidance of the two chains in a loop — 
performing the ladder sum 

The denaturation bubble under the action of the exter- 
nal force / behaves as a set of two parallel springs, each 
spring under stress //2. To find the generating function 
of the loop we form the integral 

^{zi,Z2) J dRi\dRj^G{zi,R\i,Rj^J/2)x 

G(z2,i?||,i?j.,//2). (40) 
Plugging the last line of Eq. into Eq. (gOl) yields 

T{zi,Z2)(x J d^r^ exp |-ry^l - zi/zq + yp4:X\\ 

- 2:2/20 + Z/jj/4a;|i +Z2/[|/V^} (41) 

with 

Applying this procedure to the integral in Eq. ([1T|) . we 

2;||(/) = Xii{f/2), ?/||(/) = y[|(//2), zo = zo{f/2). (42) obtain the following expression 









J ' e~'"''dw^ 







(46) 



Exchanging the orders of integration, we are left with 
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(47) 



J^(zi, Z2) cx 



In 



^1- zi/zq + y^/4x|| + y^l - Z2/Z0 + J/jj/4x|| + / 
y^l-zi/zo + yjj/4x|| + y^l - Z2/Z0 + yjj/4x|| - y\\/^ 



(48) 



Plugging zi — y/ze^^ and Z2 = \[ze into Eq. (|48|) and 
expanding in power series of d we obtain for 



y/l-\/ie'Vzo + 2/f|/4a;|| + c.c] - y,, = (49) 

5o + 2/f?/4a;|| (l-i6i- — )4-c.c. " 



where "c.c." stands for "complex conjugate." In the 
second line of Eq. we have rescaled the variable B. 

Denoting (5 = 1 — zjz^ and expanding Eq. (I49p with 
respect to 8 we end up with 



j/||/2^/xjj" + (5 - z0 + — + c.c. 



^11 



oib^d^, (50) 



where, once again, we have rescaled the variable B. 

Performing the geometrical sum in Eq. p9p and using 
the results of the expansion in Eq. ((50|) wc obtain the fol- 
lowing expression for the denaturation bubble generating 
function 



dB 



l + ylu-uln((5 + 6'2)' 



(51) 
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In Eq. (|5ip we have deformed the contour of integration 
to a unit cirlce and have assumed that the dominant con- 
tribution to the intergral comes from a small section of 
the contour between the pole of the integrand in Eq. ([?T|) 
and the branch point at 



/ //2 y^/' 

Uinh//2 



(52) 



A route to the determination of the behavior of this in- 
tegral when 6 is small is to take the derivative of the 
integral in Eq. (|5ip with respect to 6. Once the integral 
that results is evaluated, one integrates with respect to 
d to reconstruct the integral of interest. The (5-derivative 
is straightforward and results in 



Bo 



[1 + Au - u\n{5 + e-^)f 5 + 



(53) 



In the case of this integral, we can extend the upper limit 
of integration to infinity without encountering a diver- 
gence. Replacing the integration variable 6 hy y — 9/5, 
we are left with the following integral: 



dy 



1 



VS Jo (l + Au-itln(5-uln(l + y2))2 l + 



(54) 



When |uln(5| 3> 1, the dominant contribution to the in- 
tegrand in Eq. ([54|l is 



u 

71 



oo 



dy 



{\ + Au-u\n5Y ^ + y^ 

TT U 1 

'^^/6{l + Au-u\n5f 

TT U 1 



2 y/S In^ (5 



(55) 



It is now possible to integrate the expression above with 
respect to 8. This leads to the following result 



1 



Ei 



In 5 




(56) 



where Ei is an exponential integral. Expanding the above 
result with respect to b when that quantity is small, we 
find 



log(<5)^ log(<5)^ 



(57) 



which tells us that as (5 — > 0, the ladder diagram sum will 
be dominated by terms going as 



(In (5)2 



(58) 



Equation ([58)) tell us that self- avoidance leads to a log- 
arithmic modification of the contribution of the "unre- 
stricted" denaturation bubble. 



IV. THE NET EFFECT OF STRESS ON THE 
MELTING TRANSITION. 

According to the Poland-Scheraga model 0, one 
can construct the grand partition function of the dena- 
turing chain by taking the geometric sum of sequences 
of intact portions and denatured bubbles. If we call the 
grand partition function of an intact chain G(z) and the 
grand partition function of a bubble n(z), then the over- 
all grand partition function is 



G{z) + G(z)n(z)G(z) + G(z)n(z)G(z)n(z)G(z) 
_ G(z) _ 1 
~ 1-G(z)n(z) " G{z)-^-^{z\ 



(59) 



We model both double stranded (ds) and single 
stranded (ss) chains of DNA as freely jointed chains with 
different unit lengths, to take into account the greater 
flexibility of a ss chain. 

The grand partition function of the double stranded 
segment of the molecule is 



(60) 



where the summation runs over monomeric units. The 
critical fugacity of the ds chain Zds in Eq. ([50]) is given 
by 



dds/ld 



(61) 



where dds is the distance between adjacent base pairs in 
ds segment. Ids is the Kuhn length of the ds segment. 
Ag — gss — gds is the difference of Gibbs free energies per 
base pair. The second factor in Eq. (j6ip is associated with 
the configurational entropy of the ds chain (ref. Eq. (0). 

For the generating function of the denaturation bubble, 
according to Eq. we have 



n(z) 



yjl- z/Zs 



log(l - z/Zss) 



(62) 



where Zss is the critical fugacity of ss strand. From 
Eq. ^ it follows that 



/ Flss/2kBT 
\fimhFlss/2kBT 



2dss/ls 



(63) 



where, Zss is the Kuhn length of a ss strand, dss is the 
distance between adjacent base pairs in a ss segment. 

Strand separation occurs when the simple pole of the 
generating function of the DNA chain, Eq. ([59]) , coincides 
with the branch point of n(z) function, which results in 
the relation 



Zds{T,F)^Zss{T,F). 



(64) 
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Solving Eq. ([64)1 numerically, we generate the phase di- 
agram curve displayed in Fig. [S] An interesting feature 



120 



100 




120 



FIG. 8: The phase diagram: critical temperature vs the ap- 
plied force F. The following parameters were used: distance 
between adjacent base pairs of ss chain dss = 0.58nm, and of 
ds chain, dds ~ 0.34 nm; persistent length of a ss chain 0.7nm, 
of a ds chain 50nm, Ag = Ah — TAs with As = 12.5fc_g, 
Ah = TcAs, Tc = 360K JJ]. 



of the phase diagram, as was also discussed in , is the 
presence of a turning point where dT/dF changes sign. 

The relation between the strain and stress for the freely 
jointed chain of link length I is given by [12] 



AL 
17 



Nl 



coth X 



Fl 



(65) 



In the limit of the weak forces the chain's response to 
stress is linear, with effective spring constant 



cx 



(66) 



Since this spring constant in inversely proportional to the 
chain's length, a weak force aligns longer chains more eas- 
ily than shorter ones. Thus the double stranded state of 
the DNA is more favorable than the denatured one in the 
weak force limit. As the force increases in strength the 
difference of stretching per base pair of ds and ss chains 
decreases and becomes negative. This is the point where 
dT/dF changes sign. At large forces, when the molecule 
is stretched nearly to its contour length, the denatured 
state is energetically more favorable. The distance be- 
tween neighboring base pairs in an ss chain is greater 
than in the ds state due to the unstacking of base pairs. 
Breaking a base pair makes the molecule longer, thus 
reducing its potential energy, —FL. 

Next, we turn to the thermodynamic behavior of the 
system at the melting transition. In the termodynamic 
limit (TV — > oo) the free energy of the chain is dominated 
by the singularity of the generating function closest to 
the origin. Thus for the free energy per monomeric unit 
of the DNA chain we have 



Fn/N ^\ogz*{T,F), 



(67) 



where z* (T, F) is the pole of the DNA generating func- 
tion, Eq. dSH). 

To explore the behavior of the heat capacity of DNA 
in the close vicinity of the phase transition at Tc{F) we 
approximate the denominator of Eq. (j59[) , 



^l-z/z,siT,F) 



ZdsiT,F) 



log l-z/z,,(r,F) 



(68) 



by expanding the critical fugacity for ds segments of the 
chain in the first order of the reduced temperature t = 
[T — Tc)/Tc while keeping the force F as a parameter. 



Zds{T,F) ^ ZdsiT,,F) + T, 



dzds{T, F) 



dT 



t (69) 



Thus the denominator D[z) becomes 



D{z) = 1 - 



Zds{Tc,F) 



A{F)t + a 



y/1 - z/zs 



l0g(^l - z/Zss'j 



with 



A{F)^T,{F) 



dzds{T,F) 



dT 



(70) 



(71) 



To simplify even further, we drop the term 1 — z/zc in 
Eq. ([70]) in comparison with A{F)t and rescale the re- 
duced temperature 



i=A{F)t/a, 
thus obtaining for the denominator D{z) 

D(z) (xi+y/A/ log A 

where we have denoted 

A = l-z/z,,(T, F) 



(72) 



(73) 



(74) 



The pole of function D(z) can be found numerically by 
solving the equation 



A/ log A 



= -t 



(75) 



In Fig. ^ we display the behavior of A vs reduced 
temperature t. Taking the double derivative of A{t) in 
Eg. ([75)1 with respect to reduced temperature t, we find 
for the heat capacity per monomeric unit: 



(76) 



with A being the solution to Eq. ((75|) . 

As can be seen from the Fig. [HI in the vincinity of the 
critical point A cx with exponent a > 1. Therefore 
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FIG. 9: Numerical solution of Eq. (|75|l . In the vicinity of 
the critical point the DNA's fugacity approaches the critical 
fugacity as a power law It]", with exponent a > 1. 



one may conclude that the heat capacity of a DNA chain 
under applied stress behaves as 



Ccx(log|i|)^ 



(77) 



If we denote zi and Z2 fugacities associated with ds and 
ss segments in the denominator of the DNA generating 
function, 



D(zi,Z2) ^l-z,/zdsiT,F) 



y/l-Z2/Zss{T,F) 

[\ogil-Z2/zss{T,F))Y 



(78) 



then the fraction of denatured base pairs O can be de- 
termined 1231 as 



e 



d 



1-e 



D{0,z) dz 



D(0,z) 



(79) 



In Fig. [To] we plot the fraction of denatured base pairs 
for various values of the cooperativity parameter a. 



V. INFORMAL AND PRELIMINARY 
INVESTIGATION OF THE RESULTS OF 
PULLING THE CHAINS SEPARATELY DURING 
THERMAL DENATURATION 

In this section we briefly and incompletely consider the 
consequences within the Poland-Scheraga model of DNA 
melting of equal and opposite forces applied to the two 
separate chains, as illustrated in Fig. [T|3. The discussion 
below can be seen as a kind of "mean field" extension 
of the PS model to take into account the concentration 
of the potential energy associated with stretching into 
the intact base pairs. A general description of the pro- 
cess is relatively straightforward. As denaturation bub- 
bles form, the fraction of base pairs that remain bound 
diminishes, and the burden on those intact base pairs 
of maintaining the integrity of the duplexed form con- 
sequently increases. At some point, the forces tending 
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FIG. 10: Fraction of denatured base pairs O vs temperature 
T at the force F=65pN for various values of the cooperativity 
parameter a. Dotted line corresponds to a = 1, dashed line 
- to (7 = 0.5, dotted line - to ct = 0.1. Again, the physical 
parameters are the same as in Fig. [8] 



to rupture the molecular bonds between those pairs ex- 
ceed the threshold required to force the base pairs apart; 
the consequent accumulation of forces on the remaining 
bonded pairs induces a cascade of bond disruption, re- 
sulting in total strand separation. 

In the remainder of this section, we make use of the 
unmodified Poland Scheraga model to describe the sta- 
tistical mechanics of melting. We append an energy term 
associated with pulling the chains apart, assuming that 
the bonds act as Hooke's law springs under the action of 
an external force. If one pulls on a spring with spring 
constant k using a force with magnitude F, then the 
equation for the total energy in the spring is 



1 



kx^ — Fx, 



(80) 



where x satisfies the equilibrium equation 
F = kx. 



(81) 

Inserting the solution to Eq. ([8T|) into Eq. (j80|) we end up 
with the following expression for the total energy in the 
spring: —F'^ /2k. Now, start with a set of N duplexed 
base pairs. The original spring constant will be equal to 
k[)N . If N2 of the bases have been separated, the net 
spring constant is reduced to k{){N — N2). Let the force 
be equal to NFq . This results in a total energy equal to 



F^N^ 
2{N-N2)' 



(82) 



The partition function is now weighted by the factor 
exp [Fq iV^/(iV — A^2)] , where we neglect the tempera- 
ture denominator in the exponent, which can, in any case, 
be absorbed into a redefined Fq. 

As for the remainder of the partition function, we ob- 
tain it from the Poland Scheraga grand partition function 

-1 



(zq - zi+ct + d{zo - Z2Yy 



(83) 
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In the above expression, zi is the fugacity of the intact 
pairs and Z2 is the same quantity for the denatured ones. 
The first step in calculating the partition function is to 



extract the coefhcient in Eq. 
for that coefficient 



of 



{N-N2) 



{zq + ct + d{zo — Z2)P) 



Wc find 



(84) 



where we have discarded a term of order one in the ex- 
ponent. The next step is to find the coefhcient of 
in Eq. (|84p . Here we use the steepest descents method, 
looking for the extremum of 

-{N- N2) In (zo + ct + d{zo - Z2f) - N2 In Z2. (85) 



Taking the derivative with respect to Z2 of Eq. ([85|) and 
setting it equal to zero, we end up with the following 
relationship 



N2 



N- N2 



Z2dp{zn - Z2Y ^ 
zo + ci + d{z[) — Z2Y 

X{t,Z2). 



(86) 



We now introduce the variable r = N2/N. Then the 
equation above tells us that r = X{t, Z2)/(l + X{t, Z2)). 
We can use Eq. ((86)) as an equation giving us Z2 in terms 
of t and r, or as an equation for r in terms of t and Z2. 
As our final step, we plot the total free energy, given by 



{N - N2) In (zo + ct + d{zo - Z2)P) 



+ iV2lnz2 



2{N-N2) 



(87) 



as a function of A^2- We do this by first dividing ([87| by 
N, then replacing N2/N by r and, finally, plotting the 
result as a function of r. Figure [TT] shows what this pro- 
cedure leads to a range of external forces, Fq. The plots 
are for various values of reduced temperature, t. Note 
that the curves indicate an infinite minimum of energy 
at r = 1 corresponding to the lack of a limit on the lower 
bound of the potential energy when all bonds are broken 
(A'2 = N) and the strands are pulled infinitely far apart. 

The curves in Fig. [Tl]are intended to be indicative, il- 
lustrating the presence of a metastable state of relatively 
intact DNA, corresponding to a minimum with respect to 
r when the reduced temperature is sufficiently negative 
and a disappearance of the metastability as the reduced 
temperature increases. 



VI. CONCLUSIONS 

We have considered the denaturation of double 
stranded DNA under the applied stress within Poland 
and Scheraga (PS) model [2, 3]. In its original formula- 
tion, the PS model predicts a melting phase transition. 
Depending on the value of an exponent p, the phase tran- 
sition can be first or second order. In turn, the this expo- 
nent is modified by the self-avoidance of the denaturation 
loop with itself and the rest of the chain [5|- 




D.Dl 

r = N2/N 



■D.Dl 
-D.D2 
■D.D3 



FIG. 11: Free energy as a function of r for various values of 
t. The lower curves are for lower (more negative) t's. Note a 
local minimum for sufficiently negative t and an infinite drop 
in energy to the right. 



We find that external stress transforms the action of 
self-avoiding interactions. In particular, because of the 
externally applied stress, it is energetically unfavorable 
for the loop to interact with the rest of the chain, while 
for the self-intesections within the denaturation bubble, 
only a subset of interaction configurations gives a sizable 
contribution to the loop's generating function. This re- 
sults in a new analytical form for the loop's generating 
function, Eq. ([S^ . As a consequence, the phase transi- 
tion acquires a new signature; the heat capacity of the 
DNA chain behaves logarithmically in the vicinity of the 
phase transition (see Eq. (|76p ). 

Another way to apply external stretching forces is to 
pull on opposite strands of the DNA, as shown in Fig.([T]). 
In this case the applied tension energetically favors the 
completely denatured state at any temperature, leading 
either to immediate strand separation, or, under certain 
conditions, allowing the incompletely denatured DNA to 
exist as a metastable state. If there is a transition from 
this metastable state to a condition of complete denatu- 
ration, it is of the spinodal type. 

The key feature of biological DNA that complicates 
discussions of its thermal denaturation is the inherent 
inhomogeneity of its structure, the result of the fact that 
the base-pair sequence is necessarily non-uniform. As 
has been noted previously, and as we demonstrate in Ap- 
pendix [Cl the Harris criterion applies with regard to the 
relevance of this inhomogeneity, which can be treated as 
effectively random in the context of the denaturation pro- 
cess. Given the logarithmic modification of the specific 
heat at the melting transition induced by self-avoidance 
when there is melting under stress, we find that inho- 
mogeneity is relevant in three dimensions and will thus 
alter the asymptotic behavior of the system at and in 
the immediate vicinity of the transition. Precisely what 
sort of change the inhomogeneity induces has been inves- 
tigated [12, H^, [131 ; further study will no doubt prove 
valuable. 
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APPENDIX A: THE DISTRIBUTION 
FUNCTION. 

The end-to-end_probability distribution in Eq. ^ can 
be simplified to [2g| 



Pn{R, f) 



drj r/ e 



N 



.(Al) 



The integral in Eq. (jAl[) can be evaluated with the use 
of the saddle point approximation [21 ]. 



Expanding the function P/v {R, F) about the maxi- 
mum of $ to the second order in p ~ p* we arrive at 



N 



/sinh/y (N 
PNipJ)ccl^—^j exp|-.^^_^^ 



(2) (2) 

where the functions $p||P|| (/) and ^plp±{f) are 



(2) 
PllPlI 



(2) 

P±P± 



dp\ 



P — sinh^/ ' 
P 



P„,Pa 



/coth/-l 



(A9) 



(AlO) 



.(All) 



APPENDIX B: HIGHER ORDER CORRECTIONS 
TO THE STATISTICS OF AN EXCLUDED 
VOLUME POLYMER UNDER TENSION. 




PNiR, F) 



where p = R/Nl, pn = R\\/Nl is the component parallel 
to the force, and x is solution of the equation: 



coth X — — = p. 

X 

In the following we will use the notation 

Fl 



f 



(A3) 



(A4) 



To find where the function Pn{R,F) peaks we look for 
the extremum of the argument of the exponential in 
Eq. jAS]) 



HPhP±) = f p\\ +iog 



sinhi 



-px 



with 



Thus 



p = ^p\ + pi- 



d'^jpW^pi) 
dp± 

d<^{p\\:P±) 

dp» 



^ pl = 0, 



(A5) 



(A6) 



(A7) 



= 



p| =p* =coth/- -. (A8) 




'-'4 




FIG. 12: Upper left: diagram contributing to the self-energy 
of the propagator for a double strand of DNA subject to an ex- 
ternal stress. The numbers index the vertices in the diagram. 
Lower right: expanded view of the vertex labeled 4. The ar- 
rows indicate the "sense" of the strand, and the heavy dashed 
line is the effective repulsive interaction associated with the 
self-intersection occurring at the vertex. 



The essence of the analysis of higher order corrections 
and their effect on the analytical structure of the gener- 
ating function of a strand of DNA under stress is based 
on the argument that infrared singularities in the gener- 
ating function arise from the large-distance behavior of 
the real space propagator and, in particular, from inte- 
grations over large separations in the multiple integrals 
represented by perturbation-theoretical diagrams. Con- 
sider, for instance, the diagram pictured in Fig. 1121 Con- 
sider, in particular, the vertex of that diagram indicated 
by the dashed circle. The contribution to the singularity 
structure of that diagram of interest to us results from 
the integration over the position of that vertex. The in- 
fluence of the external tension on the statistics of the 
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walk is encoded in factors of the form exp[/ • (rj — rj)], 
where i and j are the indices of the vertices at the "head" 
and "tail" ends of a propagator line. As a given diagram 
consists of a single path broken into segments, each of 
which is a propagator line, the products of those factors 
is the overall exp[/ • R], where R is the displacement 
vector from the "tail" to the "head" of the DNA strand. 
This means that vertex 4 is associated with the product 
of factors 

exp [/ • (r4 - ri + ra - r4 + r4 - ra + r2 - r4)] . (Bl) 

This tells us that the force has no direct influence on the 
integration over the real space position of vertex 4. The 
principle contributions to the integration will be of the 
general form 



exp[— \/6/aZ \/l — z/zo + f^a/Ql 
|r4 - ri| + 2|r4 - rs] + - r2|)], 



(B2) 



where multiplicative factors that are independent of r4 
have been omitted. When r4 is sufficiently large, the key 

contributions go as e" 



actors that are polynomial in r^. Whatever singularity 
results from the integration over r4 arises from the coeffi- 
cient Y^l — z/ zq + f^a/Gl in the exponent, which means 
that the contribution to the singularity structure in the 
complex z-plane is at Zc in Fig. [51 which is, as noted 
previously, further from the origin than the principle sin- 
gularity in the generating function for the stress-affected 
propagator. 

All vertices to be integrated over will be of the kind 
discussed immediately above, with the exception of the 
vertex labeled 2 in Fig. [121 the analysis of which parallels 
the discussion in Section fllBI Thus, for this diagram — 
and we believe to all orders in perturbation theory — the 
effects of self-intersection are asymptotically negligible. 



APPENDIX C: DEMONSTRATION OF THE 
RELEVANCE OF DISORDER TO THE MELTING 
TRANSITION IN DNA 

The demonstration in this Appendix should be seen as 
a recapitulation and extension of the discussion by Mon- 
thus and Garel [l3] . In order to assess the effects of inho- 
mogeneity on the statistics of the melting transition, we 
make use of the fact that the "disorder" associated with 
the distribution of base pairs is quenched rather than an- 
nealed, in that the base pairs are effectively frozen into 
place and do not rearrange in response to free energy gra- 
dients. In this case, the appropriate disorder average to 
take is over the free energy, or the logarithm of the par- 
tition function, rather than the partition function itself. 
We then make use of the standard and useful identity 



In a = lim - 



1 



In order to evaluate the right hand side of (jCip we 
raise the partition function of denaturing DNA to the 
n**^ power. Then, consider the disorder at a given site, 
which is assumed to be on average equal to zero, and 
which has a gaussian distribution. Figure [13] illustrates 
the result of the disorder averaging of this disorder. The 




FIG. 13: The averaging process. The left hand side of the 
figure depicts the n lines generated by taking the n**^ power 
of the partition function of a melting strand of DNA, the large 
dots representing the disorder on a particular site. The right 
hand side of the figure depicts the result of the averaging of 
this disorder over an ensemble. 



dots on the lines indicate the disorder that is averaged. 
This disorder represents a departure from the average 
value of, say, the binding energy of the base pairs. 

To complete the process associated with disorder aver- 
aging, it is necessary to calculate the number of ways in 
which n lines can be picked out and paired. This is just 
n{n — l)/2. If we divide the result of this by n, and then 
set n equal to zero, we end up with the result depicted 
in Figure 1141 The diagonal line in the figure represents 




FIG. 14: The result of the disorder averaging of the -n}^ power 
of the grand partition function of denaturing DNA, after di- 
viding by n and setting n = 0. 



division. The figure stands for the combination 



goiz) 



2 ■ 



(C2) 



(CI) 



There are two powers of the unaveraged grand partition 
function in the denominator because the original expres- 
sion had two powers of the partition function in the pair 
that is impurity averaged, and the total power of g is, in 
the end, equal to zero. 

We now consider the mathematical structure of the ex- 
pression A in (|C2[) . It corresponds to two strands of DNA 
on each of which there is an "impurity potential" on one 
of the sites. Because the quantity being averages is the 
grand partition function, the length of the strands be- 
yond the common site containing the impurity potential 



14 



is variable. That is, for each strand in that pair, we sum 
over all lengths beyond the point at which there is an 
impurity potential. Not only that, each sum is indepen- 
dent. This means that we are left with two strands that 
have the same number of sites to the left of the impurity 
potential, but for which the number of sites to the right 
can vary. Because of this, the expression A corresponds 
figuratively to a digram like the one shown in Figure [T51 
multiplied by go(z)'^. This multiplicative factor is can- 
celled by the in the denominator of (|C2p . Thus, the 
ratio in (|C2p is pictorially represented by the diagram in 

Fig. \m 



FIG. 15: The diagrammatic result of dividing A by go{z)'^ as 
in Eq. ((02)) . 



As the next step, we turn to the expression represented 
by Fig. [151 The important characteristic of this expres- 
sion is that the number of sites in each of the two strands 
leading up to the common disordered pair is equal. How- 
ever, we are still performing a grand partition function 
sum. This means that it is necessary to work out how to 
extract the subset of terms in the two-strand sum corre- 
sponding to the same number of base pairs in each strand. 
That is, it does not suffice to simply take the product of 
two grand partition functions together. Progress can be 
made once one accepts that the principal contributions 
arise from the pole in the partition function that is closest 
to the origin. 

We start with the Poland-Scheraga propagator 0, Q 



1 



-t + {Zc - Z)P 



(C3) 



We assume that the system is below the melting tem- 
perature; the temperature is thus explicitly taken to be 
negative. The pole in the propagator in (|C3|) is at 



z = Zc-f/P 

= zc{t). (C4) 

Now, let z = Zc{t) — A. Then, the propagator becomes 

1 1 



~t + {zc ~ Zc + t^/P + A)P -t + t{l + At-i/p)p 

1 



pAt^-^/p + 0(A2) 
^(i-p)/p 



pA 



(C5) 



This tells us that the residue at the pole goes as t'^^ p)Ip , 



The proper combination of the two partition functions 
corresponding to the replicas pictured in Fig. [15] is rep- 
resented as the following sum 



f2{l-p)/p 



Zc{t) 



^2(l-p)/p 



1 



p^ Zc{tY ^ Z? 
l2{l-p)/p I 

^2 {z,{t) - z)iZcit) + Z)' 



(C6) 



Given that we are interested in the behavior of this ex- 
pression in the immediate vicinity of the singularity at 
z = Zc{t)we are left with a lowest-order contribution to 
the effect of disorder on the partition function that is 
proportional to 



^2(l-p)/p 
Zc{t) - Z 



(C7) 



If we assume that this amounts to a modification of the 
argument of the logarithm in the free energy, we have a 
new free energy that goes as 

In (z,(t)-z + A<2(i-P)/f) 
= In (^Za - t^/P + Ai2(i-P)/P (C8) 

In order that the disorder have a vanishingly small ef- 
fect on the behavior of the free energy, one demands that 
when t is very small, ^^(i-p)/? ^ i^/p qj. 2(1 — p)/p> p. 
Solving for p, we see that this is equivalent to requiring 
p < 1/2. Now, the free energy of the model is controlled 
by the behavior of Zc{t), which means it goes as t^/P. 
Thus, the specific heat — the second temperature deriva- 
tive of the free energy — goes as t~2+i/p = ^-a^ This tells 
us that the specific heat exponent is given by a — 2 — 1 /p. 
li p < 1/2, then a < 0. Thus, in order for the disorder to 
be irrelevant, we must have a < 0. Otherwise, the dis- 
order cannot be ignored. Note that this benchmark for 
the relevance of disorder to the thermodynamics of the 
melting transition is consistent with the Harris criterion 

[Til. 



The effect of logarithms 



Given the results of Section llll Al we see that the power 
p that one associates with the melting transition in the 
presence of a force is p = 1/2, with a logarithmic correc- 
tion. That is, the corresponding propagator goes as 

(C9) 



-i+(^,-z)i/2/(ln(z,-z))2- 
The pole of this propagator will be at z — Zc — S, where 



t. 



(CIO) 
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An iterative solution to the equation (jClOp yields 



This tells us that the residue goes as t{\ntY 



= 16<^(lnt)''. 



(Cll) 



We then write z = Zc — S ~ A and calculate the residue 
by expanding the following in A: 



-t + 



(16t^(lnt)4 + A)i/2 
(ln[16t2(lni)4 + A])2 

1 A 



-t + t[l 
i(ln<)4 



32 t^ilnty 



(C12) 



Once again, we process the disorder term as in (jC6|) - 
(jCSp . incorporating it into an altered argument of the 
logarithm, and we have a modified free energy going as 



ln(zc(t) - z + At^ilntf) 



(C13) 



It is clear that the contribution to the argument of the 
"disorder" term will dominate the shift in the pole in the 
ordered model when t is small. This means that disorder 
is relevant — if just barely so — for stress- modified melting 
of DNA. 
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